Spectral stability of vortices in two-dimensional 
Bose-Einstein condensates via the Evans function 
and Krein signature 

Richard KoUar 
Department of Applied Mathematics and Statistics 
Faculty of Mathematics, Physics and Informatics 
Comenius University 
Mlynska dolina, 842 48 Bratislava, Slovakia 
E-mail: koUar fmph.uniba.sk 

Robert L. Pego 
Department of Mathematical Sciences 
and Center for Nonlinear Analysis 
Carnegie Mellon University 
Wean HaU 6113, Pittsburgh, PA 15213 

Abstract 

We investigate spectral stability of vortex solutions of the Gross-Pitaevskii 
equation, a mean-field approximation for Bose-Einstein condensates (BEG) 
in an effectively two-dimensional axisymmetric harmonic trap. We study 
eigenvalues of the linearization both rigorously and through computation 
of the Evans function, a sensitive and robust technique whose use we 
justify mathematically. The absence of unstable eigenvalues is justified a 
posteriori through use of the Krein signature of purely imaginary eigenval- 
ues, which also can be used to significantly reduce computational effort. In 
particular, we prove general basic continuation results on Krein signature 
for finite systems of eigenvalues in infinite-dimensional problems. 



1 Introduction 

Background. Since the experimental creation of Bose-Einstein Condensates 
(BEC) in alkali vapors in 1995 [51 [13], BEG are one of the most active areas 
of modern condensed-matter physics. A general overview of the subject can be 
found in [121 El] i and particularly in the review book [42] . In the Hartree-Fock 
mean-field approximation, BEC are modeled by the nonlinear Schrodinger equa- 
tion (NLS) with a non-local nonlinearity. A traditional simplification, replacing 
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the non-local interaction potential with a localized short-range interaction pro- 
portional to the delta function, leads to the Gross-Pitaevskii equation 



ih^pt 



(1) 



where h is Planck's constant, M is the atomic mass of atoms in the condensate, 9 
is an azimuthal angle in cylindrical coordinates, and g is an interaction strength 
parameter. The total number of particles N in the condensate is given by the 
integral 



and is conserved during the evolution of the system. Equation ([T]) is a nonlinear 
Schrodinger equation with a cubic nonlinearity (focusing or defocusing depend- 
ing on whether the interaction is attractive or repulsive, respectively) and with 
a spatially dependent trap potential V^(x) stationary in a frame rotating with 
frequency O about the vertical axis. A rigorous mathematical justification of 
the Gross-Pitaevskii model for the BEG ground state under various conditions 
directly from many-body Schrodinger equations was done in a series of papers 



From the point of view of nonlinear waves, the interesting phenomena is that 
the Gross-Pitaevskii equation, similarly to some other nonlinear Schrodinger 
equations, supports the existence of various types of solitary wave solutions. 
In the two-dimensional setting we will study in particular, there are vortex 
solutions which have the form 



where r, 9 are polar coordinates, m is vortex degree, /i is the vortex rotation 
frequency (physically, chemical potential) and w{r) is the radial vortex profile. 

Problems of stability for vortex solutions to various forms of nonlinear Schrod- 
inger equations have drawn much attention in recent years. Related questions 
for models arising from nonlinear optics, micromagnetics and Bose-Einstein con- 
densation have been considered extensively in the mathematical and physical 
literature. For recent work concerning spectral stability questions for various 
types of matter- waves, including vortices, vortex rings, multi-poles, soliton and 
vortex necklaces in the presence of magnetic traps and optical lattices, see for 
example [36l EH EH |42l [43] . Rigorous mathematical results on these questions 
are rather few, however, due to the strong nonlinearity and complexity of the 
system. 

This work. In the present work we study the spectral stability of a sin- 
gle two-dimensional axisymmetric vortex trapped in an axisymmetric harmonic 
trap. For this simple well-studied physical setting, we develop an approach 
that involves a combination of analytical and numerical tools which allows us 
to obtain reliable results for large particle number, well into the Thomas-Fermi 
regime (e.g. N ^ 10^ atoms of ^^Na). It is important to note that in the present 
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of Lieb et al. gSJ-Ili. 



e-*^*e*"^w(r), 
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axisymmetric setting, rotation of the trap does not influence the dynamic sta- 
bility of vortices, as the rotation term can be removed by transforming to an 
appropriately rotating coordinate frame. 

Analytical results. Our study involves a number of analytical results that we 
prove in sections 5 and 6 concerning the spectrum of the operator one obtains 
by linearizing about a vortex solution. In section 5, we prove that due to the 
harmonic trapping potential the essential spectrum of this operator is empty — 
hence the spectrum consists entirely of isolated eigenvalues of finite multiplicity. 
The real part of any eigenvalue satisfies an explicit bound depending only on 
the (non-dimensionalized) frequence /i and degree m of the vortex, namely 



The eigenvalue problem breaks into an infinite system of coupled pairs of or- 
dinary differential equations (ODEs) for azimuthal Fourier modes indexed by 
j S Z. As is known from |18| , only finitely many of these are relevant for possible 
instability, namely the ones satisfying 



We construct a globally analytic Evans function [31 (TS] associated with each of 
the ODE pairs, whose zeros correspond to the eigenvalues. These results extend 
the approach of [54' for focusing-defocusing nonlinear Schrodinger equations to 
handle a harmonic trapping potential. 

One of the weaknesses in the numerical investigations in [54, was that one 
could not be confident that all unstable eigenvalues were detected, due to the 
absence of any bound on their imaginary part. No such bound is available in 
the present problem either. Nevertheless, we will show how one can indeed 
account for all possible instabilities through use of Krein signature (the sign of 
the linearized energy of associated eigenmodes) . 

The key property of Krein signature that makes it useful is its invariance 
under continuous variation of parameters such as the standing- wave frequency 
and the size of the condensate. In section 6 we extend results that are well-known 
for finite-dimensional systems to establish a general continuation property for 
any family of operators that is resolvent-continuous, a natural notion of weak 
continuity that one expects to hold in many applications to infinite-dimensional 
systems. Preservation of a definite signature is proved for any finite system of 
imaginary eigenvalues for such a family. By consequence, the only way eigenval- 
ues can leave the imaginary axis is through collision with eigenvalues of opposite 
Krein signature. 

There are three reasons why Krein signature is extremely helpful and a pow- 
erful tool in the present work. First, it allows us to explain the collisions (or 
avoided collisions) of eigenvalues found in our numerical computations. More- 
over, in combination with numerical plots, it provides a numerically convincing 
a posteriori justification that there are no unstable eigenvalues outside a certain 
fixed box in the complex plane. Finally, it can be used to significantly reduce 
the amount of necessary computation, as we shall discuss in section 7.2. 



|ReA| < 3{fi~m). 



(3) 



< IjI < 2m. 
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Numerical methods. In order to locate zeros of the Evans functions using the 
argument principle, we design and implement a numerical method somewhat 
more robust than that used in [51|. We use a path- following technique and 
multiple shooting to compute the nonlinear wave profile, and a rescaled exterior- 
product representation of the Evans function to handle problems of rapid growth 
and numerical dependence. A numerical path- following technique for radially 
symmetric profiles was also used by Edwards et al. [14]. 

Numerical results. In agreement with previous studies in the physical litera- 
ture [T51 inZl [Sni IM] ) we find a singly-quantized vortex [m — 1) spectrally stable 
while the stability of multiply-quantized vortices (with m = 2 and 3) depends 
on the diluteness of the condensate, with alternating intervals of stability and 
instability as Ng varies. Pu et al. |60) in their analysis use a different numerical 
method (finite elements) which is a priori less reliable and sensitive than the 
approach used here. Moreover, our results account for the appearance of all 
instabilities through the tracking of all eigenvalues of negative Krein signature. 
The presence of unstable eigenvalues (eigenvalues with positive real part) for 
certain parameters is also corroborated by direct simulation of time-dependent 
dynamics based on the splitting scheme proposed in [6]. 

Symmetries. In the computations one observes a special set of eigenvalues 
which remain constant under variation of standing- wave frequency and conden- 
sate size. These eigenvalues arise from the symmetries of the Gross-Pitaevskii 
equation, as we discuss in an appendix. One symmetry is particularly in- 
teresting — a breather boost, which we found first described in |59| . This is 
self-transformation of the Gross-Pitaevskii equation related to the Talanov lens 
transformation, involving a time-periodic dilation of space with an appropriate 
radial phase adjustment. This symmetry corresponds to eigenvalues ±2ia; of 
the Gross-Pitaevskii equation with the harmonic potential V'(x) = ia;2|x|2, lin- 
earized about a central vortex, with mode shapes corresponding to infinitesimal 
breathing oscillations. 

Related literature. Let us now discuss some known results on stability which 
are relevant to our problem. In general, two different concepts of stability are 
distinguished in the literature: energetic and dynamic stability. A solution 
is energetically stable if it minimizes an associated energy functional within a 
certain class of functions. 

The simplest energetic stability approach, where the minimization takes into 
account only single vortex solutions with different charges, indicates that a 
high enough trap rotation frequency can eventually stabilize a vortex of any 
degree [5]. On the other hand, without external trap rotation, the energy of a 
single multi-quantized vortex of charge m is larger than the energy of m singly- 
quantized vortices, and thus multi-quantum vortices are believed to be unstable. 
The total energy in this case also depends on the relative location of vortices as 
they tend to form regular hexagonal arrays in harmonic traps. 

A mathematical framework for a rigorous variational approach was discussed 
by Aftalion and Du Their method for effectively 2D condensates is parallel 
to the Ginzburg-Landau theory of superconductors. In [66] the authors claim 
that sufficiently fast rotation in combination with a strong pinning potential is 
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capable of making even multi-quantum vortices energetically stable. 

A detailed rigorous analysis was conducted by Seiringer [65], who studied 
regimes when a vortex solution can be a global energy minimizer. He proves that 
for any < D, < flc there exists tuq (independent of an interaction potential) 
such that all vortices with charge m > ma are energetically unstable; i.e., they 
are not global minimizers (ground states) of the energy functional (see also 
[30]). Moreover, he proves that all multi-quantized vortices, m > 2, become 
energetically unstable for a large enough value of the chemical potential of the 
condensate. Finally, he proves that symmetry breaking of the axisymmetric 
vortex solution is inevitable for any m, even for a singly-quantized vortex for a 
large enough interaction strength, since no ground state is an eigenfunction of 
the angular momentum. The symmetry breaking of a one-dimensional ground 
state is also demonstrated by a dynamical systems analysis in [33] in the case 
of a double-well trapping potential. 

Energetic stability provides a sufficient condition for dynamic stability — 
ground states of the energy functional are nonlinearly orbitally Lyapunov stable, 
i.e., if the initial data are "close" to the ground state solution then the perturbed 
solution remains "close" to the ground state solution for all times. (See [33] for 
a detailed dynamic stability study of a one-dimensional model and a sketch of 
the proof of well-posedness for the initial- value problem for the Gross-Pitaevskii 
equation.) As dynamic stability need not imply energetic stability, however, it 
may not be possible to draw conclusions on dynamic instability directly from 
considerations of the energy functional. The study of linear or spectral stability 
can be helpful since one may detect possible instabilities due to the presence of 
eigenvalues in the right half-plane. 

In the physical literature, Garci'a-Ripoll and Perez-Garci'a [18] and Pu et 
al. [BD] have studied linear stability of single multi-quantized vortices using equa- 
tions equivalent to those here. The numerical results of [60] agree substantially 
with those of the present paper. In [18] the search for instability is restricted 
only to the so-called anomalous modes - modes with a negative linearized energy. 
As pointed out in |31j for example, these modes are not intrinsically unstable in 
the sense that some dissipation mechanism must be introduced into the system 
for them to become relevant. The numerical techniques used in [TB] and [SD] 
rely on Galerkin-type approximations. A finite-temperature generalization is 
further studied in [TTj . 

Use of the Krein signature as a tool to study stability of nonlinear waves 
has appeared recently in various studies, see [9] [Ml EH EH and references 
therein. Skryabin in [57] (also see [B8j) studies a binary mixture of trapped con- 
densates using such information. Kapitula et al. 37 study stability of various 
types of matter-waves including localized vortices. Their perturbation argu- 
ment, combined with topological methods based on Krein signature, describes 
in detail transition to instability in the limit of weak atomic interactions. Fi- 
nally, Kapitula and Kevrekidis [351 ISS] have studied Bose-Einstein condensates 
in the presence of a magnetic trap and optical lattice, making efficient use of 
information on the Krein signature of relevant eigenvalues. 

Organization of this paper. First, Sections 2 and 3 introduce notation and re- 
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call background results regarding the Gross-Pitaevskii equation and vortex solu- 
tions. Section 4 contains most of our analytical results, concerning linearization, 
essential spectrum, bounds on eigenvalues, and reduction to ODEs. Moreover, 
we establish a precise asymptotic description of eigenfunctions necessary for 
construction of the Evans function. The Evans function itself is constructed in 
Section 5. In Section 6, we discuss the Krein signature. We describe in detail 
our numerical procedure and discuss the numerical results in Section 7. Finally, 
in an Appendix we discuss symmetries and boosts of the problem and relate 
them to the eigenvalues which do not change as the standing wave frequency 
varies. 

2 The Gross-Pitaevskii Equation 

The behavior of low-temperature Bose-Einstein condensates (EEC) trapped in 
the harmonic potential V{x) rotating with the angular velocity ft about the 
0-axis is well described by the time-dependent Gross-Pitaevskii equation. The 
wave function ?/;(x, satisfies ([T]) (in three dimensions) with trapping potential 
y (x) = V{x, y, z) given by 

The interaction strength parameter g is 

9 = 93D = -j^ , 

where a is the s-wave scattering length |:45 . The term fide corresponds to the 
angular momentum $7 • (r x V) of the condensate caused by a rotating frame of 
coordinates. The total number of particles in the condensate N is given by the 
integral 

/ li^fdx^'^N, (5) 

and is conserved during the evolution of the system. For disc-shaped (pancake) 
traps {ojI ^ ujI, ujI ^ u)y) it was justified [6] that the system is well approxi- 
mated by a planar two-dimensional reduced model. The equation ([T|) formally 
does not change; one only needs to set 

V{x) 
9 

where Ux = (jJtr and ujy = ujtr^tr- For purpose of numerical investigations in 
this work the same values of parameters were used as in [60 : a condensate 
consisting of atoms of ^■^Na is considered with a = 2.75 nm, lUz — 2Tr x 200 Hz, 
uJtr = 27r X 10 Hz {lUz ^ ^tr), M — 10^^^ kg and the Planck constant h = 



V2d{x) = -MujI{x^ + XW) 



92D — 93D 



/Mcjz 



V 2TTh 



1/2 
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6.6261 X 10 Js. A similar set of parameters used in the experiment as cited 
in [32 [66] with a ^^Rb condensate is: M = 3.81 x IQ-^^ kg, a = 5.77 nm, 
Lo = ujz = ojtr = 27r X 200 Hz. In these experiments the number of particles 
was approximately N = 2 x 10^, horizontal and vertical condensate sizes were 
i? = 20 /zm, L = 10 /im and the temperature Tc = 1 ^K. Both ^^Na and ^'^Rb 
represent alkali gases with a repulsive interaction potential. As an example of 
an attractive interaction, ''Li with a = —1.45 nm can serve. Note that the 
parameter a can be tuned via Feshbach resonance. 

In this work the following assumptions will be made. We assume that the 
magnetic trap is axisymmetric (Xtr — 1). Thus only two-dimensional wave 
functions of the form tp = ip(t,r,9) will be considered. The symmetry of the 
trap also eliminates dependence of stability of axisymmetric vortices on the trap 
rotation. Hence, we will set $7 = 0. Moreover, although the model includes both 
attractive and repulsive interaction interparticle potential (sign of the nonlinear 
term) , for simplicity only the more interesting case of repulsive potential will be 
considered here (some results concerning stability in transition between repulsive 
and attractive potential can be found in [60i ) . 

To nondimensionalize the equation ([l} we use the following scaling 



t=(l/wtr)t', X = y/h/MuJtrx' , ^ = y/ HuJtr / \g\ Tp' ■ 

The Gross-Pitaevskii equation is then expressed (dropping the primes) as 

iVt^-^AV^ + ir^V'+li^PV (6) 

with 

/•oo /'27r 

/ / \i^\'^rd0dr = K =\g\NM/h^ . (7) 
Jo Jo 

The energy functional is given by 



= / / 

Jo Jo 



rdO dr . (8) 



Note that the Thomas- Fermi regime (THl US] Na/do >> 1 (here do is the 
mean oscillator length, do = ^JH/MijJq and wo is the mean trap frequency Wq = 
) corresponds to if — >■ oo since K = 2\a\N y/2TTMujz/h, i.e., 



K = i:p2V2^ . (9) 

do J 

This is the limit under which Lieb and Seiringer [45] justified the Gross-Pitaevskii 
energy functional to be a good approximation for the iV-body quantum system. 
Note, that the only free parameter which stays in ©-([T]) is K , the i^-norm of 
the wave function ip. 
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3 Vortex Solutions 



In this section we describe the structure of vortex solutions [51] to the Gross- 
Pitaevskii equation ([6]), of the form 

V'(<,r,6l) =e-*'^*e™^w(r), (10) 

and describe the numerical procedure which allows us to approximate them 
with high precision. Here m is an integer, and in this paper it suffices to always 
assume that m > 1 due to reflection symmetry. 

The radial profile function w{r) of a vortex solution satisfies the equation 

1 IPTh Q I I O / \ 

— Wrr Wr H + T W + '^\w\ W — 2/iW = , T > 0. (11) 

We will require that the radial profiles be non-negative and spatially localized, 
i.e., satisfy the boundary conditions 

w{r) is bounded as r — > 0+, w{r) — > 0^ as r — > oo. (12) 

We note that the boundary condition as r — > 0+ implies w(r) — > 0+ for m > 1. 

Existence. The existence of positive solutions to (jlip for some /i, correspond- 
ing to any given K > in can be proved using a well-known variational 
argument 65 . One minimizes (jH]) among functions with the spatial dependence 
in (jlOp. with m and K fixed, and a minimizer may be found that is positive. 
For any positive solution of (jlll) with finite energy, necessarily fi > m + 1, since 
multiplying (|lip by 2T:rw and integrating in r yields 

2 I / ™ I 2 \ 2 



hK>t:J yw:^ + +r'^ j j rdr > {m + l)K. (13) 

The last inequality follows since the integral is minimized at a positive solution 
of 

1 TO^ 



— Wr 



2 w + r^w ~ 2fiw = , (14) 

which is ((TT|) linearized at zero. Analysis of this equation (see below) yields 
fl = m + 1, w = cu'q™'' from (jl8p below, where c is constant. 

The following proposition describes a global bound on any vortex solution 
with positive profile. A proof can be found in [44] (except for the statement 
that /i > m + 1). The bound ([T3l) was also proved in [65] . 

Proposition 1. Let w(r) be a finite-energy positive solution to Ul\) . satisfying 
ijJiip . where m > is an integer. Then > m + 1, and w{r) is increasing on 
(0,i?) and decreasing on (i?, oo), for some R e (m/y^S/I, y^S/i). Moreover, for 
all r > we have 

\w{r)\'^ < m. (15) 
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Figure 1: The radial vortex profile of a multi-quantized vortex, m — 2, for the 
dimensionless parameter /i « 35 corresponding to N ^ 10^ particles of ^'^Na. 
The quadratic profile in the Thomas-Fermi regime for the same parameters is 
also plotted (dashed line). Detailed (quadratic) behavior close to the origin is 
on the inset. 



The asymptotic behavior of a vortex profile can be determined directly from 
(fTTj) (see Fig. [T]) . It is not difficult to see that as r — >• 0"*" the equation has the 
same character as the linear Schrodinger equation ((T4| (one argues as in [29]). 
and 

w{r) ^ dor"^ for m > 1, 

for some positive constant dg- As r — >■ cxd, the nonlinear term for a localized 
solution becomes negligible and the linear equation (IT4l) is again a good approx- 
imation. The proof that the positive solution ■w{r) to ([TT]) approaching as 
r — >■ cx) satisfies 

w{r) = ©(r^-^e-"'/") 

is also given in detail in |44) . 

The goal of this paper is to study spectral stability of solutions to (|lip 
both analytically and numerically. Naturally, for a careful numerical stability 
investigation it is crucial to obtain very precise numerical solutions of ([TT]) first. 
The approach used here is based on path-following along a branch bifurcating 
out of the trivial solution w = and is similar to the one used in [H] . 

Bifurcation. The bifurcation (and later stability) analysis requires detailed 
information about the localized solutions to ([14]). This equation has two in- 
dependent general solutions — products of a polynomial, a decaying Gaussian 
and a confluent hypergeometric function. The exact solutions (taking jl = fJ. 
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henceforth) are 

The confluent hypergeometric functions M{a, h, x) and U{a, b, x) are, in general, 
independent solutions to xf" + {b — x)f' — af = [T]. Their asymptotics as 

r — > 0^ and as r — >■ oo is, respectively, 

w,,,{r) ^ r-{l + 0{r^)), (r) ^ i^^.--- V^^ (l + 0(l/r^)) , 



The Wronskian of W(i)(r) and ?i;(2)(^) is given by 

W^Ki)>-(2)) = -'^^^- (16) 

The only possibility for W(i)(r) (and similarly for W(^2){t)) to satisfy the bound- 
ary conditions at both ends is when the Wronskian (jl6p vanishes. This happens 
if |r((rn + 1 — /^)/2)| — oo, so m + 1 — fi = — 2n, n a nonnegative integer. 
Therefore a non-trivial solution u'„(r) to (jl4p approaching zero as r — >■ 0^ and 
as r — > oo exists if and only if /i = /i„ , where 

fi,-^ ^ m + 1 + 2n, n = 0,1,2,.... (17) 

For /J, given by ([T7| both solutions (r) and 10(2) (r) reduce to a constant 
multiple of the single solution 

wi"^\r)=r"^e-^'/'Li"^\r'), (18) 

where L^T'^r) is the generalized Laguerre polynomial, with n the number of 
zeros of L^\r) for r > 0. The positive solution (the ground state of the 
associated energy functional) corresponds to n = 0, /j,o = m -I- 1 and w — 

The numerical algorithm designed to find solutions to ([TT|) is based on the 
following observation. It is reasonable to expect that introduction of the non- 
linear term leads to the existence of a solution branch (/i(s), Wf^{s)) bifurcating 
from the trivial solution w = for fi = fiQ. To justify such a behavior one can 
use the Crandall-Rabinowitz theorem [TTJ [53] , to prove the following Theorem 
(for details see [H]). 

Theorem 1. The solutions (/i,w) to ill]) near (/io,0), iiq — m + 1, form a 
curve 

s ^ {p.{s),w{s)) = {fiQ + t(s), sw^™' + sz{s)) , 
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Figure 2: Left panel: (log^Q N vs. /x) plot of branches of vortex solutions for 
m = 1, 2, 3 (from left to right) emerging at /xq = m + 1 for ^'^Na data. Dashed 
curve represents the number of particles in the Thomas-Fermi parabolic regime 
WTF{r) ^ \J [i- r'^/l. 

Right panel: The difference d — J I^'tfP — |^'|^ of the number of particles of 
vortex solutions and the number of particles of wtf for m = 1,2,3. 
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where s i~> (t(s),z(s)) G R x span |u'q™''| is continuously differentiable near 
s^O, t(0) = t'(0) = 0, z(0) = and 

w{s) eX = {w: e*™^w(r) G V.} , 



Note that in addition to a branch of positive solutions {p,, w) bifurcating 
from the trivial solution at /x = /io = m + 1 in a direction (0,Wq'"^) there are 
sign-changing branches of solutions (^(s), wi™^ (s)), n — 1,2,..., bifurcating 
from the trivial solution at /i = fin = m + 1 + 2n in a direction (0, wi™-*). The 
proof is analogous. 

Numerics. Hence it is possible to numerically trace the solution curve (/i, w) 
from the branching point {fio,0), see Fig. [2ja). The typical radial profile for 
m = 2 is on Fig. [T] Note that the Crandall-Rabinowitz bifurcation theory also 
provides some linear stability information for solution branches given in Theo- 
rem [T] [TTl [53] , but one must take the stability results with caution. Although 
the theory predicts stability for the bifurcating branch, it is only with respect 
to radially symmetric perturbations in (jl4p . This is not sufficient to determine 
stability of vortex solutions with respect to the full dynamics in ([6]). 

The behavior of norms relative to the norm of the parabolic Thomas-Fermi 
regime approximation is illustrated in Fig. [Hb). The first point on the ap- 
proximate solution curve is set to be {fj,Q , ewl^'^ ) , e — 0.1, which is an O(e^) 
approximation of the exact solution. Then an implementation of a predictor- 
corrector path-following algorithm [4] is used to get solutions for large values of 
the parameter /i. 

Since our stability study will require evaluations of the profile at any given 
point within the computational domain, the precision of calculations is improved 
by optimizing the already calculated profile for any given /i by a multiple shoot- 
ing procedure [69] . This allows one to achieve high precision in evaluating w{r) 
by simple integration from a nearby mesh point. Note that the calculation is 
almost independent of the size of the parameter fi since the size of the com- 
putational domain, and so the number of necessary nodes, grows very slowly. 
Therefore it is possible to reach large values of fi. Also note, that with the grow- 
ing parameter ^, the L^-norm ^ of profiles grows (Fig. ^ and hence states far 
in the physically interesting Thomas-Fermi regime for a wide range of /x's are 
obtained for a small computational cost. On the other hand, as pointed in [14] . 
for computation for a single value of /i this method has significant overhead. 

4 Linearization and Reduction to Ordinary Dif- 
ferential Equations 

The goal of this section is to derive and study the linearization of (fTTj) around the 
solutions constructed in the previous section — localized vortex profiles. The 
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linearized equations have the same form as the so-called Bogoliubov equations 
[Mj [58] commonly used in physics literature. Note, however, that the rela- 
tion between the derivation here and the physical derivation of the Bogoliubov 
equations is by no means straightforward. 

A small general perturbation of the vortex solution ^{t, r, 9) = e'''~''*"'"™*'w(r), 
where w{r) = Wn{r) for a fixed parameter /z, has the form 

uit,r,e) = e-'"' {e'"''^wir)+ev{t,r,9)) . 

Neglecting nonlinear terms in ^ yields 

ivt ^ -^A-y - fiv+ ^r'^v + 2|wpu + \w\'^ e^'^^^v . (19) 



The complex character of the equation ([191) complicates the analysis. Therefore, 
we decompose the complex wave function v as 



$i\ _ /Re V 
$2/ ^ \Im V 



The equation ([r9[) is then equivalent to the real system 



= A$ = ./ 
where 



$ , (20) 



j=(^J oj ^^''^ ^=(o -1 

To understand the dynamic stability of the vortex solution -0, we will study 
the spectrum of the operator A as an unbounded operator on L^(M^,R^) with 
an appropriate domain D{A). As is rather well-known, spectral stability of A 
(meaning absence of spectrum in the right half of the complex plane) , need not 
necessarily imply linear stability (in the sense that the zero solution of ([20[) is 
stable), nor nonlinear stability of vortices. In this paper, we avoid these subtle 
issues and confine ourselves to studying spectral stability. After all, the presence 
or absence of eigenvalues with positive real part is interesting in itself. 

The precise definition of the operator A is somewhat involved and requires 
the concept of a quadratic form 62 . Write (only formally for now) 

i^c^-^A + ir^/, L^=2\w\''I+\w\^e^^''R-^i, (21) 

where / is the 2x2 identity matrix. Then 

A^-J{L, + L^). (22) 

Then define a quadratic map 

qLc ■■ D{q) = {H^{R^, ]R2) n L2(r2^ jj2. ^2)^ ^ ^ 
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by 



1, 



The quadratic form q^c is semibounded: 



1*1 



r dr dO . 



is dense in L^(I 



Here 



Note that the space D{q) = 7J1(R2,M2) ^ i2 

since the Schwartz space is a subset of and is dense in L^(R^,R^ 
L^(M^, R^; /(r)) represents the space of functions 4' : R^ ^ R^ with the bounded 
norm 



1*1 



L2(R2^R2.^(r)) 



27r 



\^{r)\'^f{r)rdr d0 . 



The quadratic form q^c is closed if it has a closed graph, i.e., if D(q) is 
complete under the graph norm |l*||+i 



/gLc(*,*) + 11*11^2. This is true 
since both and L'^(r'^) can be obtained from by completing the space of 
functions under the and L'^{P) norms respectively. Also observe that 



C 11*1 



1*1 



L2(r2) 



>9Lc(*,*) 



I*lli2 >c 



1*1 



//I 



1*1 



L2(r2) 



for some C > c > 0. (The lower bound follows from the semiboudedness; 
the proof of the upper bound is analogous.) Then Theorem VIII. 15, pp. 278 
of [62] yields that q^c is the quadratic form of a unique self- adjoint operator 
Lc- The domain of the operator Lc, denoted by D{Lc), is dense in L^. Clearly 

H^(M.\M.^)nL^n 



l'^;r^) C D{Lc), and we have 



L,* 



1 



-A 



1 



-r^I * 



for * G D{Lc) in the sense of distributions. 

It is easy to see that the operator L^, in (PT|) is bounded on L^(R^,R^) 
Therefore the operator A ^ — J{Lc + Ly^) with the domain D{A) = D{Lc) is 
closed and densely defined in L^(R^,R^). 



4.1 Essential spectrum 

We investigate the spectrum (y{A) of the operator A given by ([2^ . regarded as an 
unbounded operator on complexified space D{A) C L^(R^,C^). The spectrum 
of such an operator in general consists of two parts: isolated eigenvalues of 
finite multiplicity form the discrete spectrum crdisc(^), and the remaining part 
— the essential spectrum iJcas{A). The latter is empty as stated in the next 
proposition. 

Proposition 2. For any m and ^ the spectrum of A consists entirely of isolated 
eigenvalues of finite multiplicity. I.e., the essential spectrum of the operator A 
is empty. 
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Proof. The proof has three steps. First, it is easy to see that the essential 
spectrum of Lc is empty by Theorem XIII. 16 on p. 120 of Reed and Simon [63] . 
Then we prove the same for JLc- Finally, the generalized Weyl theorem for 
non-self adjoint operators yields the same for JL. 

Let us prove that the essential spectrum of JLc is empty. Since Lc is a non- 
negative operator, is not an eigenvalue and therefore Lc has bounded inverse. 
Moreover, Lc only has discrete spectrum and its eigenvalues are isolated with 
the only possible accumulation point oo. Then by Theorem XIII. 64 on p. 245 
of [53] the operator L~^ is compact. Now consider the following identity: 

\I - JLc = {XL-^J-^ - I)JLc. (23) 

If A crdisc(>/ic), then j ^ ad\sc{Lc^ J^^) and the right-hand side of (|23)) has 
bounded inverse given by 

(A/ - JLc)'^ = Lc^J-\XL-^J-^ ~ I)-^ . 

Here L'^ is compact, J^^ = —J is a bounded operator, XLj^J^'^ — / is a 
compact perturbation of identity, a Fredholm operator. Moreover, since j ^ 
o'disciL^^ J^^), the operator XLc^J^^ — I has empty kernel and is invertible 
and bounded. Therefore {XI — JLc)~^ is compact. This implies that XI — JLc 
is invertible with a compact inverse if A is not an eigenvalue of JLc- 

It remains to prove that the eigenvalues of JLc are isolated and of finite 
multiplicity, which prohibits discrete spectra to be embedded in the essential 
spectrum. To show that, consider the resolvent equation 

{IX - JLc)u = f . 

which is equivalent to 

-(/ - T{X))u = (AL-ij-i -I)u^ L-^J-^f . 

The operator XLc^J^^ — / is a (multiple of) compact perturbation of identity 
and therefore it is also Fredholm. Also, it is analytic everywhere except for the 
discrete spectra of Lj^ J~^. By a general result of Gohberg and Krein [22], p. 21. 
or Kato [JT], p. 370, the set of values for which /— T(A) is not invertible is at most 
countable with their only possible accumulation point infinity. Therefore the 
eigenvalues of JLc are isolated. Also, the spectral projection on the eigenspace 
associated with a particular eigenvalue of JLc has finite dimensional range, since 
it is given by an integral 

Px^^ f{XI-JLc)-'dX 

of a compact operator. Hence the essential spectrum of JLc is empty and consist 
of isolated eigenvalues of finite multiplicity with only accumulation point infinity, 
i.e., 

acss{JLc) = . 
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Finally, we use the generalization of Weyl's theorem to non-self-adjoint op- 
erators to prove 

(TossiJLc) — acss{J{Lc + L^)) . 

It is enough to observe that J{Lc + L^) is a relatively compact perturbation of 
JLc, i.e., that JL^{\I — JLc)^^ is compact whenever A ^ ddisdJLc)- This is 
true since JL^, is bounded and {XI — JLc)^^ is compact. □ 

Let us point out that although the stability of an axisymmetric vortex in an 
axisymmetric trap does not depend on the trap rotation frequency Q, it has an 
influence on the linearized operator A. The definition of the operator and the 
proof of emptiness of the essential spectra can be adjusted to account for the 
rotation as long as |ri| < uitr- Beyond this threshold it is not clear how to define 
the operator and whether its essential spectrum stays empty in this parameter 
regime. This threshold may be significant if the axial symmetry is broken. A 
further discussion on other features in this regime can be found in Chapter 12 
[42] . On the other hand, it is easy to see that the eigenvalues (the discrete 
spectrum) of the linearized problem suffer only a shift by a purely imaginary 
number depending on rotation, so stability of this part of the spectrum of A is 
unaffected by rotation of the trap. 

4.2 Eigenvalues 

By the result above, the investigation of the spectrum of the operator A is 
reduced to study of the eigenvalue equation 

A$ = A$ , (24) 



$ + |w|2e'™^^i?$ = 0. (25) 

Similarly as in [54] it is useful to represent solutions of ([25]) in the basis of the 
eigenvectors of the matrix J: 

Since D{A) C L'^{R'^,C'^), then $± G L'^{R'^,C). Consequently, ^ has the 
form 

(^tX+^A-K^ +ti-2\wf^<^>+~\wfe''"^'^^- = 0, (26) 

(^-iA+iA-ir2+/i-2|w;|2^$_-|u.pe-2™''$+ = 0. (27) 

Using the information on asymptotic decay and a simple bootstrap argument 
one can deduce that $± g H^^^ for each fc > 0, so $± € C°° (R^ ,C) n L^{M.^ ,C). 



which can be rewritten as 



A J 



-A- 



- - 2|w|' 
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Furthermore, decompose ^±{r,0) at each fixed r mto Fourier modes (with 
shifted indices for notational ease) 



$±(r, 



E 



9 (j±m) I \ 

y± {r) • 



(28) 



After the introduction of the Fourier modes the equation (1251) transforms to an 
infinite-dimensional system of linear equations. 



The system decouples to coupled pairs for nodes y+ — y/ 



(j+m) 



,y- 



y 



iX + -Ar 



-iX + -A,, - 



(j + m)^ 1 2 ,9 



U - ™)^ 

2r2 



1 

r~ 

2 



-^-2|it;|' 



y+ 

y- 



(29) 

w\^y+. (30) 



By the symbol we denote the radial Laplace operator A^ — ~^ r'§p- 

The proper boundary conditions for this system must be determined only 
from the system itself and the property $± S L^(M^,C). Here asymptotic 
behavior of solutions to ([29|) ~(p0 | described in Theorem [2] of the next Section 
can be used. It implies that the appropriate boundary conditions are 



lim y±{r) exists and 



lim y±{r) — . 



(31) 



Therefore the eigenproblem for A is decomposed into countable many prob- 
lems 

Ljy = iXRy, y = {y+,y- f , (32) 



where 



and 







Li 



(j ± m) 
2r2 



2 1 
1 2 



2^ 



/i. 



Similarly as before the bootstrap argument gives y G C°° D ( 
associated inner product is given by 



(33) 

-,C2;r-). The 



iy,z)^ / (y+{r)z+{r) +y^{r)z^{r)'jrdr . 



On the other hand, one can argue (as in [54]) that any solution (A, y, j, m, /i) to 
(021) defines a solution (A, $, m, fi) to ^ with $ £ L'^{R'^,C^). 

Note that = L^j, a fact connected to the Hamiltonian symmetries of A. 
The next proposition analogous to !54^ summarizes the results of this section. 
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Proposition 3. A complex number A is an eigenvalue of A if and only if for 
some integer j the system of equations iSS^) have a nontrivial solution satisfying 
i31\) . An eigenfunction of A associated with an eigenvalue A has the form 

If form a solution to I13S\) for some pair (A,j), then {y-,y^) form a 

solution for (—A, — j) and {y-, y+) form a solution for (A, — j). 

4.3 Bounds on unstable eigenvalues 

At a first glance it may seem impossible to solve infinitely many systems of the 
form (j3D) . Fortunately, similarly as in [M] it is possible to restrict the index j for 
which an unstable eigenvalue may occur. The proof of the following Proposition 
can be found in [T5] but for the sake of completeness and clarity of exposition 
we provide it here as well. 

Proposition 4. All the possible unstable eigenvalues and eigenf unctions of the 
operator A must correspond to bounded solutions of \31]) -^i ^] for j satisfying 

j^O and \j\<2m. (34) 



Proof. First we prove that \j\ < 2m. The main idea of the proof of this part is 
the same as in |18) but we present it here for clarity. Assume the contrary, i.e. 
IjI > 2m, namely, |j — m\ > m and \j + m\ > to. The operator Lj is analogous 
to the operator Lc in Section 2] and thus the integrals used here are well defined. 
Then for any y ^ by integration by parts 



{Ljy,y) = 



j + m 



\y+\ 



Ifj-m 
2 



+ y+y-) - ^ {^ry+y+ + ^ry-y- 



> 



1 

2^2" 



1 



\y- 



Here we also used the inequality 

'{y+y- + y+y-)< \y+\'' + \y-? ■ 

Hence {Ljy,y) > Ew{y+) + E^{y^) where 



\drf? 



TO 



-r^ + 2\w\- 



1 2 

-r 
2 



r dr 
1 



2\w[ 



M (I2/+I 



{\dry+\' + \dry-\^) 



(35) 



r dr . 



2/i l/l^ 



r dr . 



(36) 
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The function w{r) is a non- negative minimizer of the hnearized energy (|36p 
within the family of functions /(r) with (j){r,6) — e™^/(r) G D{Lj). Since 
E^{w) = it foUows that E^{f) > for ah e*™^/ G D{Lj). Hence 

(Lj-y, y) > . 

If (A, y) are an eigenvalue and eigenfunction vector as in (j32p then also 

{iXRy,y) =iX{Ry,y) > 0. 

Therefore lA must be real if |j| > 2m. 

This result can be also interpreted in terms of Krein signature (see Section[5]): 
the signature of all eigenvalues A for \j\ > 2to is positive. 

Next, we prove the stronger result that \j\ < 2m and j ^ 0. First, let us 
consider the case j = 0. In that case one obtains directly 

{Ljy,y) > E^{y+) + E^{y^) . 

The equality is valid only if there is equality in (1351) . i.e., if y+ = — y-. Therefore 
{Ljy,y) = only if j/+ — —y- almost everywhere and Ew{y+) = 0. Since 
/(r) — w{r) is a non-negative solution of the differential equation associated 
with dMl) 

-frr - -/,. + ^/ + r^f + 2|uf / - 2/i/ - , 

the second condition holds only if ?/+ = aw where a is a complex number, and 
w — w{r) is a real ground state of the energy E^lf). Thus 

y = w{a, —a) . 

Then the eigenvalue problem ((^^ - ([50)) reduces to iXw = 0. Hence A = and 
there are no unstable eigenvalues for j = 0. Similarly, one can rule out the case 
j = ±2m. □ 

This result can be also interpreted in terms of Krein signature (see Section|6]): 
the signature of all eigenvalues A for \j\ > 2m is positive. 

One can also prove the following estimate which restricts the possible un- 
stable eigenvalues to lie in a vertical strip. 

Proposition 5. The real part of every eigenvalue of the operator A is bounded: 
|ReA| < 3niax|u;(r)p <3(p-m) <oo. (37) 



Proof. First, recall the simple bootstrap argument justifying that $± of (l26l) - 
([27)1 is C°°(M^ M^) n L2(M^ M2)_ rpj^g^^ gpjj^^ i-^e operator A to a vortex-profile- 
dependent part A^^ and independent part A^: A = Ac + A^, where 



Ac 

An,! 



J 



1 2 

2' 



-J [2\w{r)\^ + \w{r)\^e 



2 2m0J 



R] 
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Since J and R are constant matrices (bounded by 1 in the matrix norm) the 
estimate (jisp imphes that the norm of the profile dependent part is bounded 
above: 

\\A^\\l2 = \\~j{2\w\^ + |u.pe2"^"'i?)||^, < 3MH < cx) , (38) 
where II -11^2 denotes the operator norm in (R^, C^) and M(u)) — max |u'(r)p. 

re(0,oo) 

Multiply ((25|) by the smooth complex conjugate $ and integrate over to 
obtain 

A||$f = / / Ac<^-'^rdrde+ j / A^$-$rdrd6'. (39) 
Jo Jo Ja Jo 

The second term on the right hand side of can be estimated using (P5| 

/ / A^$-$rdrd6i < 3M(w)||$|p. (40) 
Jo Jo 

Finally, the real part of J^^ Ac^ ■ <i> r drdO vanishes, which can be checked 
by a simple but long integration by parts (omitted here) . The statement of the 
proposition then immediately follows by Proposition [TJ □ 

Propositions S] and [5] restrict the search for unstable eigenvalues to a finite 
number of equations (with indices < j < 2m) and to a vertical strip. Since 
one can expect infinitely many stable eigenvalues in both directions on the 
imaginary axis it is hopeless to prove that the imaginary part of an eigenvalue 
is bounded. Nevertheless, the possible number of unstable eigenvalues is limited 
(see Section in]). 

5 Evans Function 

While the finite element and the Galerkin approximation methods provide a 
fast and simple way to find eigenvalues of the problem (15^ , the Evans function 
technique [3l [15] method has proved to be the most reliable and robust in certain 
cases. This approach will be implemented here. It is parallel to [54], where the 
reader can find many details of the procedure. 

The main idea of this approach is to identify eigenvalues of the operator Lj 
as zeros of an analytic function Ej{\). First, write the system (|^ - (|5n|) as a 
4x4 system of first order ordinary differential equations: 

y'^B[r,j,\)y (41) 

where 

B = B^+B.,„, y=[y^^"'\r),dry^+^\rly^-'-\r),dry^l-^\r)f (42) 
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and 





1 







k+ 


-l/r 

















1 


Vo 





k- 





Bo 



The coefficients fc+ and fc are given by 

(i + 



r 










4 





2 

















V2 





4 


0/ 



- 2^ - 2iA . 
-2^ + 2iA. 



k+{r, A) 
fc-(r, A) 



Tlie asymptotic behavior of solutions to (j4T|) is described in the next theorem. 

Theorem 2. Jbr any A G C and m > 0, fi real, j integer, there exist solu- 
tions yf^\r) and yi°°\r), i — 1,2,3,4, to the system j[ j with the following 
asymptotic behavior, 



yooi{r) 



as ?■ — > 0+ , 
as r ~> CO . 



Here yoi and yooi, i = lj2,3,4, are independent solutions of the asymptotic 
systems 

Vo = Bo{r,j, X)yo , yoo = B^oirJ, X)yao , 



where 







1 





-i 


I -ly 





1 



and 



l+{r 



[j + mf 



l-ir) 



u - 



The asymptotic behavior of these solutions as r 
where = M + ct- = M ^ ^'^'^ as r 0+ &?/ 



oo is given by 



yooi - /2r"+ (l/r, 1, 0, 0)' , yoo2 - e"- /^r"- (0, 0, l/r, 1)' , 

g-rV2^-a+ (1/^^ o)T ^ _ e-rV2^-a- (q, 0, l/r, -if , 



yol -rlJ+"l (l,|j+m|/r,0,0)' 



2/02 ^rlJ-"! (0,0,l,|j-m|/r)' 



2/03 ^ (1, -|j + m|/r, 0, 0)^ , yo4 ^ r-^-"! (0, 0, 1, -\j - m\/ry 
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The full proof of the theorem which relies on the asymptotic theory of Coppcl 
[10] can be fomid in [44] . 

The asymptotic analysis reveals that (|^T|) has two exponentially growing so- 
lutions asymptotically equivalent to y'^^ (r) and yj"'' (^) foi' f and two expo- 
nentially decreasing solutions asymptotically equivalent to y^°°\r) and y^\r) 
for r ^ 1. Note that these particular solutions are not in any way unique. 
From now on the notation y[^\ y2^\ 2/3°°'' ^^d 2/4°°'' will always refer to solu- 
tions with the given asymptotics. The two-dimensional growing and decaying 
subspaces non-trivially intersect only if A is an eigenvalue. Their intersection can 

be detected by vanishing of the Wronskian W{r,X) = det(2/j°\ 2/2°^ y3°°^ ^l""'')- 
By Abel's formula, this determinant satisfies a differential equation W'{r) = 
Tr (B) W{r). It is convenient to remove the dependence on r by setting 

E,{X) = -r'det{y^°\r),y!f\r),yir\r),yi"^\r)) . (43) 

This Evans function is then independent of r. 

It is evident that Ej (A) = is a necessary and sufficient condition for the 
existence of an eigenvalue. A different representation of the Evans function is 
more convenient for analytic and numerical use, however, to avoid problems such 
as maintaining linear independence of solutions for extreme values of parameters 
and at mode collisions. It is also desirable to ensure global analyticity of the 
Evans function so that the presence and location of eigenvalues may be studied 
by means of contour integrals via the argument principle. 

An alternative way to construct and evaluate the Evans function involves 
introducing the adjoint system 

z' = ~zB{r,j,X). (44) 

The fundamental matrices Y(z) (its columns are yi) and Z(z) (with rows Zi) of 
systems (PT|) and (|33]) are related by ZY = I. Therefore Theorem[2](by a simple 
direct calculation of the inverse matrix) also guarantees existence of four inde- 
pendent solutions of (|44|) as z[^\ Z2°°\ ^3°°^ ^.nd z^^ such that the matrices 
with columns and with columns yl°°^ satisfy = /. 

One can also easily deduce the asymptotic behavior of z'^°°'^ as r ^ 00. Fur- 
thermore, a simple calculation [54] shows that 



Constructing Ej{\) in this way still involves maintaining the independence of 

particular solutions yf^ and zf^\ however. A quite simple idea to overcome this 
difficulty has been used in a number of earlier works including |54| . Instead of 
considering the system (j4ip one can construct a larger 6x6 system for exterior 
products of solutions to (|4ip . This exterior system has a unique solution of 



maximum growth rate given by ^^2'* = v'l^ ^ Vt'^ ^^d a unique solution of 
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maximal decay rate given by y^^' — yl°°' A . Corresponding statements 
hold for the adjoint system. Evaluation of the Evans function is then given by 
the simple formula 

E,iX) = ziT^-y[lK (46) 

It is important to realize that the Evans function given by (|46)) is analytic in 
C and it is solution- and spatially- independent. On the other hand, the values 
are the same as given by (|43|) and (j45|) . 

These Evans functions have the following symmetries. The proof is the same 
as in [51] . 

Proposition 6. For all integers j and complex numbers A G C, 

Particularly, Ej{X) is real for X purely imaginary and EqIX) = Eq{—X) = Eq{X) . 



6 Krein Signature 

In this section we establish some general basic continuation results on Krein 
signature for finite systems of eigenvalues in infinite-dimensional problems, and 
identify the eigenvalues of negative signature for the linearized Gross-Pitaevskii 
equation 

A typical context [24l |39] in which Krein signature arises is in study of 
linearized Hamiltonian systems 

vt = JLv (47) 

on a Hilbert space X with inner product (•,•), where v £ X, L: X ^ X is 
symmetric, {Lu,v) = {u,Lv), and J : X ^ X is invertible and skew-symmetric, 
{Ju,v) = —{u, Jv). We are interested in cases when L is unbounded and not a 
positive operator and has a finite number of negative eigenvalues. In particular, 
we will apply the results of this section to the operators L — Lj, J = J = —iR, 
see p2|) . Note that individual operators Lj that we consider here do not have 
all the symmetries of the Hamiltonian operator L = Lc + of (|2ip . 

To define Krein signature of a discrete eigenvalue A of JL (or more generally, 
a finite collection of such eigenvalues), consider the restriction of the energy form 
(L-, •) to the associated generalized eigenspace. If this restriction is positive or 
negative definite, the Krein signature of the eigenvalue is positive or negative, 
respectively. In the case when the restricted energy is indefinite, the Krein 
signature of the eigenvalue is indefinite as well. Note that it is easy to see that 
the Krein signature of each eigenvalue off the imaginary axis is zero, since if 
JLu = Xu with A 7^ —A, then (Lu, u) = since 

A( J^^M, u) = {Lu, u) = (u, Lu) = (u, A J"^u) = -A( J"^u, u). (48) 
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Let us point out that in [3D1[72] (also see [331 [SD]) it is rigorously proved for 
finite-dimensional systems that the only way eigenvalues of a system depending 
on a parameter can leave the imaginary axis (as a quadruplet, since the symme- 
try of the problem forces two complex conjugate pairs of the purely imaginary 
eigenvalues to collide simultaneously) is via a collision of two purely imaginary 
eigenvalues of opposite signature. This behavior is generic - the passing of 
two eigenvalues of opposite Krein signature on the imaginary axis is an event 
of codimension one (R. Kollar & P. Miller 2010, unpublished), as a particular 
quantity involving eigenvectors of the associated eigenspaces must vanish. 

For simplicity, in what follows we will always assume J is invertible. Then 
Xu = JLu is equivalent to Lu = XJ~^u, and thus 

{Lu,u) ^ X{J-\,u) . (49) 

For this reason it is sometimes more convenient to use the real quadratic form 
(iJ~^u,u) instead of {Lu,u). 

6.1 Finite systems of eigenvalues 

In this subsection we formulate precise statements of some key properties of 
Krein signature that apply to unbounded operators of the type we consider. 
A key fact that makes the Krein signature so useful for continuation problems 
is that for finite systems of imaginary eigenvalues, it can never be (positive 
or negative) semi- definite without being definite. This is a consequence of the 
following. 

Lemma 1. Let X he a Hilbert space, and let L be a symmetric and J a skew- 
adjoint operator on X with a bounded inverse J"^. Let S C iM &e a finite set 
of discrete purely imaginary eigenvalues of JL and let Xi be the corresponding 
spectral subspace (the span of all generalized eigenvectors for eigenvalues in Sj. 
Then the quadratic form (ij~^u,u) is non-degenerate on Xi. 

Proof. This result is well-known in the finite-dimensional case (see 72, p. 180]) 
and the proof here is not very different, based on the spectral decomposition 

X^Xi®X2. 

of the underlying Hilbert space into the finite-dimensional Ji-invariant subspace 
corresponding to E and its Ji-invariant spectral complement. The spectrum of 
JL\xi is S, and the spectrum of JL\x2 contains no point of E. (See Theorem 
III-6.17of gl].) 

We claim that {J^^u,v) — whenever u e Xi and v £ X2. This is true 
because, if u is a generalized eigenvector for some A e S, then (A — JL)^u = 
for some n, and if w e X2, then w = {X— JL)~"w G X2, and one checks that 

{J~^u, v) = ( J^^M, (A - JLYw) = ( J"^(A + JL)"u, w) = 0. 

Now, if the form {iJ^^u,u) is degenerate on Xi, it means there exists a non- 
trivial u S Xi such that {J^^u,v) = for all v £ Xi, hence for all v £ X. 
Hence J^^u = 0, sou = 0,a contradiction. □ 
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A simple corollary for simple eigenvalues immediately follows. 

Corollary 1. Let X , L, J he as above. Let A G iM 6e a simple isolated eigen- 
value of JL with corresponding eigenvector u. Then {ij^^u,u) is non-zero, and 
if Xy^O, so is the Krein signature sgn{Lu,u) . 

In the application we will consider, the operator L depends continuously (in 
a sense we will make precise) on a real parameter /i. By the corollary above, 
the only possibility for a continuously varying eigenvalue A = A(/i) of JL{fi) to 
change its Krein signature is for it to collide with another eigenvalue, or to cross 
zero. As a simple consequence of (|49|. a simple eigenvalue crossing zero will flip 
its Krein signature to the opposite. (If the change is from negative to positive 
this may correspond to symmetry-breaking instability [50 .) 

Corollary 2. Let X be a Hilbert space, and let J be a skew- symmetric operator 
on X with bounded inverse. Suppose s i— >■ L{s) is a family of symmetric oper- 
ators on X, for s in an interval I C M. Assume that for s G I, {X{s),u{s)) is 
a continuous family of isolated simple purely imaginary eigenvalues and corre- 
sponding eigenvectors for the problem 

X{s)u{s) — JL{s)u{s) . 

Then iX{s){L{s)u{s), u{s)) has the same sign for all s with A(s) 7^ 0. 

Proof. The continuous quantity {iJ^^u{s), u{s)) is real and does not vanish for 
any s e X (by Corollary [Ij . Then 

iX{s){L{s)u{s), u{s)) = X{s)^iiJ-\is), u(s)) for all sel. 

The result follows since A(s)^ < 0. □ 

The previous two results guarantee that the Krein signature of a purely 
imaginary eigenvalue does not change unless the eigenvalue crosses the origin 
or collides with other eigenvalues. Also, Theorem [5] implies that the operators 
JL = iRLj we will consider have only discrete eigenvalues of finite multiplicity, 
so only finitely many eigenvalues may collide at one point, for any value of the 
parameter fj,. 

Such colliding eigenvalues will form a finite system of eigenvalues in the sense 
of Kato (see section IV-3.5 of [41 ), with a corresponding spectral projection 
X I— > Xi{^) that varies continuously with /i. If no eigenvalue of negative or 
indefinite signature is present in this family before collision, then the signature 
is positive definite and must remain so at collision and after. Then all eigenvalues 
must remain on the imaginary axis after collision — no pair of eigenvalues can 
bifurcate off the axis, because the signature of such a pair is indefinite. 

This is well known for finite-dimensional systems [5D]. To justify these state- 
ments about continuation of Krein signature for unbounded operators, we should 
first define an appropriate notion of continuity. The following definition is essen- 
tially related to the concept of convergence of closed operators in the generalized 
sense of H], see Theorem 2.25, Section IV-2.6. 
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Definition 1. Let s T{s) be a family of closed operators in X , defined for s 
in an interval I C M. We say the family T{s) is resolvent-continuous on X if 
for each sq G X, there exists ^ € C such that the resolvent map s i-> ~T{s))^^ 
is defined and norm- continuous for s in a neighborhood of sq. 

For such a resolvent-continuous family T(s), the resolvent map is actually 
(locally) continuous in s for each point ^ of the resolvent set of T(s). 

Suppose now that Sq is a finite set of discrete eigenvalues of T(so). This set 
may be continued continuously to comprise a finite system of eigenvalues I](s) 
defined on a maximal interval of existence I' C I in a standard way: Let T be 
any smooth contour (a collection of small circles, say) that contains Eq inside, 
and no other point of the spectrum of T{sq). The spectral projection 



is well-defined and continuous in a neighborhood of sq, its range Xi{s) is 
T(s)-invariant and has constant dimension, and the eigenvalues of the finite- 
dimensional map T{s)\xi{s) comprise a finite system of eigenvalues that vary 
continuously with s. The projection P{s) is independent of the choice of F, 
and consequently it can be defined continuously in this way for all s in some 
maximal interval X' C X. Evidently the maximal interval must be open. (At 
an endpoint of X' , in general one may have collisions with other eigenvalues or 
continuous spectrum, or other pathologies.) 

Theorem 3. Let X be a Hilbert space and let J be a skew- symmetric opera- 
tor on X with bounded inverse. Suppose s i-> L{s) is a family of symmetric 
operators for s in an interval X d M. such that the family JL{s) is resolvent- 
continuous. Suppose S(s) is a finite system of discrete eigenvalues of JL(s) as 
described above, with associated generalized eigenspace Xi{s), defined for s in 
some maximal interval X' . 

Then, if the quadratic form {ij~^u,u) is (positive or negative) definite on 
Xi{sq) for some sq G X' , it has the same definiteness on Xi(s) for all s G X' . 
By consequence, all eigenvalues in S(s) are purely imaginary for all s X' . 

Proof. The proof is again a rather straightforward extension of well-known re- 
sults from finite dimensions [5D1|72]. Whenever the form {iJ~^u,u) is definite 
on Xi{s) we must have S(s) C as follows from Considering the positive 
definite case, let p{s) be the infimum of {iJ^^u, u) over the unit sphere in Xi{s). 
Then p{sq) > and p is continuous on X'. If definiteness fails at some point of 
X', then p{s) > on some maximal strict subinterval of I' and p{si) = at an 
endpoint si G X' . But then E(si) C iM by continuity of the eigenvalues and so 
p(si) > by Lemma [TJ a contradiction. □ 

Note that in our application the resolvent-continuity condition on the family 
of operators is satisfied, as the domain Y C X oi the operators i(/i) is fixed, 
and the operators JL{ii) vary continuously in L{Y,X). 
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6.2 Negative-signature eigenvalues for the GP equation 



Here we identify all imaginary eigenvalues with negative Krein signature for 
11 = m + 1, when the linearized problem reduces to the case of dU]), the Gross- 
Pitaevskii equation linearized about a trivial solution. The linearized system 
P2p for given j and m then reduces to 

- ^il)u^ iXu , (ff^-™ - ^/) w = -iXv , 

where 

1 1 

" 2 2r2 + 2 ' 

see According to Proposition U we may assume j,m > 0. Since, by 

the discussion following the eigenvalues of H'' are exactly |fc| + 1 + 2n 

for n = 0, 1, 2, ... , with eigenfunctions wi'^'*, the system eigenvectors (u, v)'^ = 
{wn~^"^\o)'^ and (0, correspond to eigenvalues 

j + m + l + 2n = m + l + iX, |j - m| + 1 + 2n = to + 1 - iA . (50) 

The eigenvalues of the linearized Gross-Pitaevskii equation for /i = m + 1 then 
come in two families for j > 0: 

A = -.(, + 2.), ,= ^(-^- + 2n) (0<,<TO), ^^^^ 

I «(— 2to + j + 2n) (to < j) , 

for n = 0, 1, 2, . . .. It is easy to determine the Krein signature of these eigenvalues 
since in the different cases we find 



{Liu, 




= 


— iJ.I)u, u) ~ 


iX{u, u) = 


{L{u, 


vf,iu,vf) 




- ^^I)v,v) = 


—iX{v, v) 


{L{u, 


vf,iu,vf) 




- f^I)v,v) = 


— zA(ti, v) 



respectively. Note that according to Proposition!?] we can restrict our search for 
possible unstable eigenvalues to the modes with < j < 2m. In particular, this 
means that when /i = to + 1 , the only eigenvalues with negative Krein signature 
are 

for TO = 1: A = — i (j = 1), 

for TO = 2: X = -i {j = l), -2i{j = 2), -i (j = 3). 



7 Numerical Methods and Results 

In this section we first describe in detail the way we evaluate the Evans function 
and determine the presence and location of eigenvalues. We then show how the 
results on Krein signature from Section [6] are used to justify the finding of all 
unstable eigenvalues and interpret eigenvalue collisions. Based on our analytical 
and numerical results we also suggest an improved numerical method to detect 
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Figure 3: Contour F used for counting of eigenvalues. Due to the symmetry of 
the Evans function, the computation was restricted to the thick contour in the 
right-half plane. 



unstable eigenvalues by tracking eigenvalues of negative Krein signature. De- 
spite we have not used the method in the present work we believe that such an 
algorithm offers a significant reduction of computational cost in a wide variety 
of numerical studies of stability problems that use the Evans function. In the 
further subsections, numerical results are discussed separately for singly- and 
multi-quantized vortices since the stability diagrams (diagrams of stable and 
unstable eigenvalues) reveal different patterns. 

7.1 Evans function evaluation 

To attain high precision and stability for all computations, exterior products 
were used throughout for numerical evaluation of the Evans function. Note that 
exterior products may be avoided, however, if one uses the algorithm introduced 

in [mill]. 

As easily seen from the asymptotic descriptions in Theorem [2l the behavior 
of solutions of the system (|^T|) is significantly different for r ^ 1 and r ^ 1. 
Consequently it is useful to rescale the solution during the integration process 
over the interval (0,oo). (The implementation approximates this interval with 
[£:,i?], where £ = 10~^ and R = R{fi) is set in such a way that the vortex 
solution is negligible at r = i?. This is chosen as an increasing function 

with i?(0) — 5, i?(35) — 25.) The aim is to rescale in such a way that the 
matrix i?(r. A, j) (and the solution) remains appropriately bounded throughout 
the integration. The details of the rescaling used here can be found in |44j . 

The presence of unstable eigenvalues is detected by contour integration using 
the argument principle, similarly as in [54], for j's restricted by Proposition |4l 
The algorithm adaptively calculates the argument of the Evans function Ej{X) 
along an approximately rectangular contour F which encloses a bounded region 
in M^. The region is pictured on Fig. [3] Note that Proposition [5] allows us 
to confine the integration to the right half plane (the total argument change 
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is twice as large) and also reduces the set of j's for which the calculation is 
necessary to positive values. Moreover, Proposition [5] restricts the location of 
unstable eigenvalues to a vertical strip. Naturally, it is not possible to perform 
numerical calculations without imposing also some vertical bound for the region 
enclosed by the contour. The vertical bound in the implementation is chosen is 
such a way that the behavior of the stable eigenvalues becomes predictable. We 
justified that all the unstable eigenvalues are included a posteriori by examining 
the Krein signature of eigenvalues. 

Stable eigenvalues on the imaginary axis are, thanks to the symmetry of the 
Evans function, zeros of the real- valued function E(X). Hence one can plot that 
real function and determine the location and multiplicity of its zeros within a 
finite interval. In the actual implementation this is done automatically — one 
first interpolates the real function by a cubic spline and then uses the Newton 
method for locating zeros. In a small neighborhood of a possible double zero, a 
very fine mesh was used to resolve any ambiguity. 

The total number of eigenvalues of Ej{X) enclosed in the region is then 
determined by the difference Usu = '^s+u — Jt-s between the total number of 
eigenvalues inside F, 

1 f E' (X^ 1 f 

=2^^i ^""^ ^^^^ 

and the number of (stable) eigenvalues on the imaginary axis, including their 
algebraic multiplicity. If risu is not equal to zero it must by Proposition |6] be 
an even positive integer and corresponds to twice the number of pairs of stable 
and unstable eigenvalues (A, A) enclosed within F. 

The precise location of unstable eigenvalues can be theoretically determined 
by the generalized argument principle. Higher moments 

« = i:Af = 5i-/^A'||^.A (52) 

give the sum of fc-th powers of positions of all eigenvalues enclosed by F. Given 
the approximate location of eigenvalues on imaginary axis and number of eigen- 
values off the axis, the problem reduces to solving a set of nonlinear equations. 
This is particularly efficient in the 2, where only si is necessary. 

Unfortunately, even in this case, the numerical error involved can be signifi- 
cant, with a major contribution coming from a finite difference approximation 
of Ej{X). Therefore the obtained values are considered only approximate. The 
calculated location is further used to construct a smaller contour Tg which lies 
solely in the right-half plane and encloses only a single eigenvalue. The presence 
of a zero of Ej{X) inside a smaller contour is again justified by the argument 
principle. Its location is then determined by the generalized argument principle. 
This process can be repeated a few times until a desired precision is attained. 
In the implementation the threshold for a precision was set up to be 10~^. 

Note that this method of locating eigenvalues does not allow one to calculate 
the eigenfunction directly; for that another method must be used. 
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7.2 Numerical uses of Krein signature 

Here we show how the Krein signature serves to provide evidence that no unsta- 
ble eigenvalues are missed in the numerical computations. Moreover we suggest 
a very efficient numerical algorithm based on the theory in Section 6. Finally, we 
discuss an approach that does not require homotopy and avoids path-following 
completely. 

In our numerical approach we consider the linearized eigenvalue problem 
for a large range of the parameter /x, continuously connecting the linearized 
eigenvalue problem about a vortex solution to a linearized eigenvalue problem 
about a trivial solution. (See [H] for a similar approach). The latter problem has 
only purely imaginary eigenvalues, and their Krein signatures were determined 
in Section |6^ We plot the location of all imaginary eigenvalues in a bounded 
interval on the imaginary axis which contains the continuation of all negative- 
signature eigenvalues for the whole range of /i considered (/i £ (m + 1, 35)). By 
calculating the Evans function on the contour T and on an appropriate portion 
of the imaginary axis, we indirectly track signature changes using the theory of 
Section 6, and can thus infer restrictions on possible departures of eigenvalues 
from the imaginary axis. This tracking process provides strong evidence that 
there are no unstable eigenvalues other than the ones we find. 

For a significant reduction in computational effort, we propose a new numeri- 
cal method that avoids the computation of contour integrals for Evans functions. 
To identify all unstable eigenvalues, it is only necessary to perform the following 
calculations: 

• evaluate the Evans function on the imaginary axis close to zero to detect 
any eigenvalues crossing zero; 

• follow the eigenvalue branches starting at negative Krein signature eigen- 
values at the reference value of the continuation parameter (in our case 
start at fi ~ m + 1 and follow a total of 2m — 1 different branches), to 
detect any collision with positive Krein signature eigenvalues; 

• follow any branches of eigenvalues that bifurcate into the complex plane. 

We emphasize that the theory presented in Section 6 justifies that no unstable 
eigenvalues can be left out. No contour integration is necessary at all since one 
can just use a path- following algorithm (a boundary- value solver). Moreover, 
global analyticity of the Evans function is not needed, since only zeroes of the 
real Evans function on the imaginary axis need to be found initially, and one 
may use a continuation algorithm [27[ 161] to track them. 

A negative feature of the homotopy technique is that it has a large overhead 
if one is interested in only a small set of values of the parameter. Therefore 
we examine ways of avoiding this overhead by directly evaluating the Krein sig- 
nature of a given eigenvalue on the imaginary axis. This is not feasible with 
the approach of evaluating Evans functions using exterior products, since the 
eigenfunction is unavailable and it must be calculated separately by a different 
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method. On the other hand, the restriction of the Evans function to the imagi- 
nary axis is a real function. This means that only a continuous representation 
of the Evans function is required and problems of rapid growth and numerical 
dependence are much less severe. By consequence, it is not necessary to use 
exterior products to evaluate the Evans function, and it can be evaluated ei- 
ther (i) by calculating a relevant Wronskian of solutions to the linear eigenvalue 
equations which makes it easy to compute the eigenfunction, or (ii) one can 
use a simple continuation algorithm proposed by Sandstede [17]. By either of 
these methods, it is possible to evaluate the Krein signature directly, which then 
allows one to determine the number of unstable eigenvalues. 

Then for an eigenvalue problem of the form JLu = Xu the following method 
may be used. First, determine (somehow, by using oscillation theory, for ex- 
ample) the number ritotai of negative eigenvalues of L. If the number is zero, 
there are no unstable eigenvalues of JL. If ritotai > 0, the number gives the 
total number Uu of pairs of unstable eigenvalues of JL plus the number Ug of 
stable eigenvalues with negative Krein signature [311 301 [53] (also see [IS]). If 
the present method (the Evans function) detects a number of pairs of eigenval- 
ues off the imaginary axis n„ = ritotai, there are no other unstable eigenvalues. 
If < ritotai , perform a calculation of the Krein signature for eigenvalues on 
the imaginary axis. If + = ritotai there are no other unstable eigenvalues, 
otherwise, increase the area of search and repeat the whole process. 

7.3 Singly-quantized vortices m = 1 

Similarily as in [1 ^ 1371 IBOl 168] we found singly-quantized vortex, m = 1, to be 
spectrally stable for all values of the parameter /i investigated, fi ^ {m + 1, 35), 
corresponding to number of particles N G (0, 10^) for the data for ^^Na given 
in Section [2] By Propositions H] [3] and [6] it suffices to study only j = 1. For 
illustrative purposes the location of the stable eigenvalues (/i vs. Im A) for 
j ~ 0, 1, 2, 3, is plotted in Fig. 2] For the sake of clarity only eigenvalues with 
I ImA| < 5 are shown. 

For small values of /i, close to ij^q ^ m + 1, the eigenvalues are close to the 
eigenvalues of the reduced uncoupled linear problem neglecting the \w\'^ depen- 
dence. As fi increases, certain eigenvalues remain constant: a double eigenvalue 
A = and simple eigenvalues A = ±2i for j = and simple eigenvalues A = ±i 
for j = 1. These eigenvalues originate in the symmetries and boosts of the 
Gross-Pitaevskii equation and are present for every m (see the Appendix). 

7.4 Multi-quantized vortices m > 2 

The eigenvalue diagrams for multi-quantized vortices with m — 2 show more 
complexity. The radial vortex profile for /i w 35 is illustrated in Fig. [TJ 

In the case m = 2 the possible unstable eigenvalues may appear for |j| = 
1,2,3. The modes j = and j = 1 demonstrate the same features as in 
the case of the singly-quantized vortex with the same constant eigenvalues: a 
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Figure 4: Stable eigenvalues of the linearization about a singly-quantized (to = 
1) vortex solution, /i e (to + 1,30): the eigenvalues corresponding to modes 
j = 0, 1 (left panel) and j ~ 2,3 (right'^anel) . 



double eigenvalue A = and simple eigenvalues A = ±2? for j = and simple 
eigenvalues X = ±i for j — I (see Fig. [5] (a)). 

A different behavior appears for modes j — 2 and j — 3, as shown on 
Fig. [S] (b) and Fig. [71 For j = 3 there are no unstable eigenvalues present and 
the stable eigenvalues do not collide but rather diverge from each other when 
they approach each other. Collisions do occur for j ~ 2 and cause instability. 
A collision of two stable purely imaginary eigenvalues produces a pair of stable 
and unstable complex eigenvalues symmetric with respect to the imaginary axis. 
After an increase of /x these two eigenvalues return to the imaginary axis and 
split to two purely imaginary eigenvalues as illustrated in Fig. [51 This "collision 
- split - collision" process ("bubbles of instability" [SOI) repeats regularly for 
the whole range of fi studied. 

This behavior is caused by the presence of a single eigenvalue with negative 
Krein signature [SUl HHj- The imaginary part of this eigenvalue is decreasing 
with increasing /x, and on its way it encounters eigenvalues with the opposite 
signature. After each collision the eigenvalues split off the imaginary axis and 
become eigenvalue pairs with zero Krein signature symmetric relative to the 
imaginary axis. Reversibility of this process suggests that the eigenvalues come 
back to the imaginary axis and the process repeats itself (for a larger parame- 
ter /i). This indicates a surprising thing: transitions to instability for larger /i 
may happen at a large frequency (Im A large) , and therefore there is no hope to 
confine the imaginary parts of unstable eigenvalues to a finite interval indepen- 
dent of II. The behavior of the eigenvalues demonstrates strong agreement with 
[60] . This is also consistent with the results of Seiringer [65] where he proved 
that for any m > 1 for large enough fi the vortex becomes energetically unstable 
in the sense that it cannot be a global minimizer of the energy and is subject to 
symmetry breaking. Also note, that the maximum real part of unstable eigen- 
values grows very slowly with growing parameter fi. The maximum of the real 
part is much smaller than the bound we were able to obtain in Proposition [5] 

In the case j = 3 the eigenvalues have the same Krein signature and therefore 
they cannot split off the imaginary axis The eigenvalue which originates at 
A = — « for /i = 3 has negative signature at first, but after it crosses zero for /x 
close to 5, it changes its signature to positive according to Corollary [H Then 
all eigenvalues have the same signature, making splitting impossible. Instead, 
eigenvalues repel each other upon approach, as expected by [50\ . 

7.5 Further results 

The presence of exponential instability was also checked by direct simulations. 
First, an approximate eigenvector was obtained by Galerkin approximation [18] 
and then the Strang-splitting scheme 6 was used for time evolution. The initial 
perturbation of a vortex solution by an approximated eigenvector showed a good 
agreement with the expected exponential growth. 

For the case m — 3, we performed similar computations but omit details for 
brevity. There are more eigenvalues of negative Krein signature, and overlapping 
bubbles of instability. See [SD] for the analog of Fig. [SI in this case. 
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Ao = -4i 


-iA(/i) = -1.41 - 2.69Af-"-96 






Ao = 


-iX{n) = 1.38 - 7.96 A*"^-^^ 






Ao = 2i 


-iA(^) = 3.11 - 27.66/1-1-59 




,7=3 


Ao = —Si 


-iX{n) = -1.74- 5.14 






Ao = -i 


-iA(/i) = 1.73- 3.97 /i-i "2 






Ao = i 


-a(/i) =3.61-5.30/i-°-9^ 



Table 1: Approximate asymptotic behavior of eigenvalues for m ~ 1,2. 

Finally, it would be plausible to describe the asymptotics of the eigenvalues 
as /i — cx) when the condensate approaches the Thomas-Fermi regime. We 
were only able to study the asymptotic behavior of purely imaginary eigenvalues 
numerically by plotting a loglog plot of the first order differences of X{[i). We 
observed a clear linear trend implying an algebraic approach to a limit. Our 
numerical results agree with the recent analytical results in [57] for stability of 
vortex solutions in the limiting Thomas-Fermi regime (see also |17 | for stability 
of ground states). 

The approximate behavior of eigenvalues with a small imaginary part for 
m — I, j = 1,2, and m — 2, j — 2,3, is presented in Table [TJ For most 
eigenvalues, the asymptotic behavior as /i 3> 1 is well approximated by 

A(/i)=5--. 

M 

On the other hand, certain eigenvalues show different rate of convergence clearly 
distinct from {—^~^), but we do not have any explanation of this phenomena. 
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Appendix: Symmetries and eigenvalues 



The symmetries of the Gross-Pitaevskii equation and its hnearization imply the 
presence of a special set of eigenvalues. 
Phase symmetries. For any m, 

A = (for j = 0) 

is a constant double eigenvalue for all fi > fiQ. Its multiplicity comes from the 
symmetry of the Gross-Pitaevskii equation under a phase change (this generates 
an eigenvector) and under a change of standing-wave frequency (this generates 
a generalized eigenvector). A detailed discussion on these two symmetries and 
their implications for spectra is given in [54_. 

GGV boost. Due to the presence of the harmonic potential, the other usual 
invariants of the nonlinear Schrodinger equation — spatial translations — do not 
apply. Similarly, one cannot perform the typical Galilean boost. Instead, one 
can apply a boost for quadratic potentials that was described by Garcfa-RipoU 
et al. PJj. This particular boost transforms any solution ^(r,t) (r ~ {x,y)) to 
a new solution of the form 

*fl(r, t) = ^{r - R{t), t) exp(z0(r, t)) , (53) 

where R{t) is the path of a classical particle moving in the potential well. For 
the harmonic potential V{r) = ^r^, this requires simply Ru = — i?, so that each 
component of R{t) can be any linear combination of cosi and sini. According 
to [19], d is given up to a constant by 0(r, t) = r ■ Rt. 

As discussed in [54], any one-parameter family Ur of solutions to ([6]) gives 
rise to a solution of the corresponding equation linearized about mq, given by 

In the case of the GGV boost ([53|l . setting R{t) = r(cos t, sini)-^ gives 6{r, 9, t) = 
Trcos{9 — t). Then Ur is given by ([53]) with ^(r. t) — ■i/j{t, r, 9) as given by (jlOp 
with w{r) satisfying (|lip . Hence 

/ cos t 

u — \ . 
\smt 



, r, 0) -f- ir cos(9 — t)'ip(t, r, 9) , 
which can be also written as 



TflVO 

w'{r) CQs{9 — t) + i sin(i ~ 9) + irw cos(i 

r 



The vector ^ = (Reu,lmu)^ then satisfies the linearized equations (j^H)) and 
has the form 

/cos m0\ /cos 6'\ /cost 
V sin m9 J \ sin 9 J \ sin t 



- sin jn9 
cos m9 



w'(r) 

f cos9\ f—smt\mw{r) /cos9\ /cost^ 



\siii9 J \ cost J r \sm9j \smt 



rw{r) 
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In the decoupling to Fourier modes this sohition yields a solution to the system 
(|29l)~(l30l) for j = 1, A = i and 



y+ 



, mw(r) 
w (r) h rw(r) 



y- 



, mw(r) 
w [r) H rw[r) 



The similar solution can be also obtained for A = 



1 



, mw{r) ■ 
w [rj rw[r) 



1 



, mw(r) 
w (rj H h rw[r) 



Therefore the GGV boost implies the existence of the eigenvalues (for any m) 

X = ±i (forj = l). 

Breather boost. Finally, we describe the source of the presence of the eigen- 
values 

A = ±2i (for j = 0), 

which appear for every m > 1 and /i > /io- These eigenvalues corresponds to 
a "breathing" symmetry of the Gross-Pitaevskii equation in 2 + 1 dimensions, 
as explained by Pitaevskii and Rosch in [59^. Here we provide a somewhat dif- 
ferent perspective, showing that the symmetry corresponds to the Talanov lens 
transformation via an exact transformation to the cubic Schrodinger equation. 
Rybin et al. [64j noted that the radially symmetric Gross-Pitaevskii equation in 
2-1-1 dimensions transforms exactly to the cubic Schrodinger equation without 
potential, according to a transformation described originally by Niederer [52] 
for the linear Schrodinger equation. More generally, they noted that this trans- 
formation works in any dimension as long as the nonlinearity is critical, of the 
form [mI^/^m in n space dimensions. This symmetry was used by Carles [7] to 
study various mathematical aspects of the Gross-Pitaevskii equation. 

By transforming from GP to cubic Schrodinger, using a Talanov lens trans- 
formation [70], and transforming back, one can get a self-transformation of GP 
(see also |21|'). The symmetry that we will describe is somewhat more general, 
and works as follows. Suppose that v{x,t) is any solution of the equation 



idtv 



1 



-Av 



1 



Lu^\x\'^v- X\v\Pv = 



2 

where t e M and a; G M". Let 

u(t,x) =ae-'^l^l'/2i;(c2;,T) 

where a, 5, c and r are functions of time t. Then we find that 



1 



1 



idtu H — Au v \x\ u — "f\u\Pu = 0, 



(54) 



(55) 



(56) 



provided that b and c satisfy 



b'^b^+LU^C^ 



be. 



(57) 
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and 



a = c 



c\ 7 = \c'~''^'\ (58) 

Of course, & = 0, c = 1, always works with — oj'^ . 

The symmetry is most interesting when the nonhnearity is critical, meaning 
p — A/n. Then 7 = A. In 2+1 dimensions this corresponds to the cubic nonlin- 
earity. For ui^ = = we recover the well-known Talanov lens transformation 
for the critical Schrodinger equation: 



1-bot 1-bot 1 - bot 

With uj'^ — and constant i/^ > we get the transformation from the critical 
Schrodinger equation to GP as described by Rybin et al. [64 and Carles [7]: 

1 tanut , 

b = ly ta.n.i/t, c= , t— . (60) 

cos lyt V 

In the other direction with constant > and i^^ = we get 

uP't 1 arctanwt . , 

b— TT-^, c = , , T = . (61) 

The consequences of this transformation for the Gross-Pitaevskii equation 
in 2+1 dimensions are striking. From any solution of (j54l) . one gets from (|55p a 
solution obtained by a breather boost — a time-periodic dilation of space with an 
appropriate radial phase adjustment. See |59j for more discussion and a relation 
to representations of the group S0(2,l). 

A short calculation shows that this transformation corresponds to the eigen- 
values ±2i of the normalized linearized equations The relation to 
these eigenvalues can be seen also from a simple consideration, that this self- 
transformation produces small oscillations at exactly twice the trap frequency 
for classical oscillations in the harmonic potential iw^jxp. The exact formu- 
lae for this breathing boost also show that breathing oscillations need not be 
small. Also note that a breather boost can be applied to any solution, not only 
standing waves. 
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Figure 5: Stable eigenvalues of the linearization about a multi-quantized vortex 
solution, m = 2, fi G (m + 1, 35): j = 0, 1 (left panel) and j = 3 (right panel), 
with a detail of an avoided eigenvalue collision on the inset. 
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Figure 6: Stable and unstable eigenvalues of the linearization about a multi- 
quantized vortex solution, to = 2, /i e (to + 1, 35), corresponding to the mode 
J = 2. 
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Figure 7: Imaginary part of stable (light) and unstable (dark) eigenvalues of the 
linearization about a multi-quantized vortex solution, m = 2, /Lt € (m + 1,35), 
corresponding to mode j = 2. A detail of a bubble of instability on the inset. 
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